/* Copyright (c) Colorado School of Mines, 2011.*/
/* All rights reserved.                       */

/* SUDATUMK2DR: $Revision: 1.8 $ ; $Date: 2011/11/16 17:45:18 $	*/

#include "su.h"
#include "segy.h"
#include "header.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
"SUDATUMK2DR - Kirchhoff datuming of receivers for 2D prestack data	",
"		(shot gathers are the input)				",
" 									",
"    sudatumk2dr  infile=  outfile=  [parameters] 			",
"									",
" Required parameters:							",
" infile=stdin		file for input seismic traces			",
" outfile=stdout	file for common offset migration output  	",
" ttfile=		file for input traveltime tables		",
"   The following 9 parameters describe traveltime tables:		",
" fzt= 			first depth sample in traveltime table		",
" nzt= 			number of depth samples in traveltime table	",
" dzt=			depth interval in traveltime table		",
" fxt=			first lateral sample in traveltime table	",
" nxt=			number of lateral samples in traveltime table	",
" dxt=			lateral interval in traveltime table		",
" fs= 			x-coordinate of first source			",
" ns= 			number of sources				",
" ds= 			x-coordinate increment of sources		",
"									",
" fxi=                  x-coordinate of the first surface location      ",
" dxi=                  horizontal spacing on surface                   ",
" nxi=                  number of input surface locations               ",
" sgn=			Sign of the datuming process (up=-1 or down=1)  ",
"                                                                       ",
" Optional Parameters:							",
" dt= or from header (dt) 	time sampling interval of input data	",
" ft= or from header (ft) 	first time sample of input data		",
" surf=\"0,0;99999,0\"  The first surface defined the recording surface ",
" surf=\"0,0;99999,0\"  and the second one, the new datum.              ",
"                       \"x1,z1;x2,z2;x3,z3;...\"                       ",
" fzo=fzt		z-coordinate of first point in output trace 	",
" dzo=0.2*dzt		vertical spacing of output trace 		",
" nzo=5*(nzt-1)+1 	number of points in output trace		",	
" fxso=fxt		x-coordinate of first shot	 		",
" dxso=0.5*dxt		shot horizontal spacing		 		",
" nxso=2*(nxt-1)+1  	number of shots 				",
" fxgo=fxt		x-coordinate of first receiver			",
" dxgo=0.5*dxt		receiver horizontal spacing			",
" nxgo=nxso		number of receivers per shot			",
" fmax=0.25/dt		frequency-highcut for input traces		",
" offmax=99999		maximum absolute offset allowed in migration 	",
" aperx=nxt*dxt/2  	migration lateral aperature 			",
" angmax=60		migration angle aperature from vertical 	",
" v0=1500(m/s)		reference velocity value at surface		",
" dvz=0.0  		reference velocity vertical gradient		",
" antiali=1             Antialiase filter (no-filter = 0)               ",
" jpfile=stderr		job print file name 				",
" mtr=100  		print verbal information at every mtr traces	",
" ntr=100000		maximum number of input traces to be migrated	",
"									",
" verbose=0		silent, =1 chatty				",
"									",
" Notes:								",
" 1. Traveltime tables were generated by program rayt2d (or other ones)	",
"    on relatively coarse grids, with dimension ns*nxt*nzt. In the	",
"    datuming process, traveltimes are interpolated into shot/gephone 	",
"    positions and output grids.					",
" 2. Input traces must be SU format and organized in common rec. gathers",
" 3. If the offset value of an input trace is not in the offset array 	",
"    of output, the nearest one in the array is chosen. 		",
" 4. Amplitudes are computed using the reference velocity profile, v(z),",
"    specified by the parameters v0= and dvz=.				",
" 5. Input traces must specify source and receiver positions via the header",
"    fields tr.sx and tr.gx. Offset is computed automatically.		",
"									",
NULL};
/*
 * Author:  Trino Salinas, 05/01/96,  Colorado School of Mines
 *
 * This code is based on sukzmig2d.c written by Zhenyue Liu, 03/01/95.
 * Subroutines from Dave Hale's modeling library were adapted in
 * this code to define topography using cubic splines.
 *
 * This code implements a Kirchhoff extraplolation operator that allows to
 * transfer data from one reference surface to another.  The formula used in
 * this application is a far field approximation of the Berryhill's original
 * formula (Berryhill, 1979).  This equation is the result of a stationary
 * phase analysis to get an analog asymptotic expansion for the two-and-one
 * half dimensional extrapolation formula (Bleistein, 1984).
 *
 * The extrapolation formula permits the downward continuation of upgoing
 * waves  and  upward  continuation  of  downgoing waves.  For upward conti-
 * nuation of upgoing waves and downward continuation of downgoing waves,
 * the conjugate transpose of the equation is used (Bevc, 1993).
 *
 * References :
 *
 * Berryhill, J.R., 1979, Wave equation datuming: Geophysics,
 *   44, 1329--1344.
 *
 * _______________, 1984, Wave equation datuming before stack
 *   (short note) : Geophysics, 49, 2064--2067.
 *
 * Bevc, D., 1993, Data parallel wave equation datuming with
 *   irregular acquisition topography :  63rd Ann. Internat.
 *   Mtg., SEG, Expanded Abstracts, 197--200.
 *
 * Bleistein, N., 1984, Mathematical methods for wave phenomena,
 *   Academic Press Inc. (Harcourt Brace Jovanovich Publishers),
 *   New York.
 *
**************** end self doc ***********************************/
/* TYPEDEFS */
typedef struct SurfaceSegmentStruct {
        float x;        /* x coordinate of segment midpoint */
        float z;        /* z coordinate of segment midpoint */
        float s;        /* x component of unit-normal-vector */
        float c;        /* z component of unit-normal-vector */
} SurfaceSegment;
typedef struct SurfaceStruct {
        int ns;               /* number of surface segments */
        float ds;             /* segment length */
        SurfaceSegment *ss;   /* array[ns] of surface segments */
} Surface;

/* Prototype */
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,float **tb,
	     float **t,float x0);
  void interpx(int nxt,float fxt,float dxt,int nx,float fx,float dx,int nzt,
	     float **tt,float **t);
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t);
  void timeb(int nr,int nz,float dr,float dz,float fz,float sz,float a,
	     float v0, float **t,float **p,float **sig,float **ang);
  void dat2d(float *trace,int nt,float ft,float dt,float sx,float gx,
	     float **dat,float aperx,int nx,float fx,float dx,float nz,
  	     float fz,float dz,int mtmax,float xm,float fmax,
             int nxi,float fxi,float dxi,float angmax,float **tb,float **pb,
             float **angb,int nr,float **tsum,int nzt,
             float fzt,float dzt,int nxt,float fxt,float dxt,int antiali,
             int sgn,float **szif,float nangl,float **sigb, int verbose);
  void decodeSurfaces(int sgn, int *nrPtr, int **nxzPtr, float ***xPtr,
	     float ***zPtr);
  int decodeSurface(char *string, int sgn, int *nxzPtr, float **xPtr,
	     float **zPtr);
  void makesurf(float dsmax, int nr, int *nu, float **xu, float **zu,
             Surface **r);
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
             Surface *srf,float **szif,float *sz,float *nangl);

/* segy trace */
segy tr, tro;

int
main (int argc, char **argv)
{
	int nt;

	int nzt;
	int nxt;

	int nzo;
	int nxso;
	int nxgo;

	int ns;
	int nr;

	int is;
	int io;
	int ixo;
	int it;

	int antiali;

	int ntr;
	int jtr;
	int ktr;
	int mtr;

        int mtmax;
	int sgn;

	int nsrf;
	int *nxzsrf;

	int nxi;
	int *ng;

	off_t nseek;

	float ft;
	float fzt;
	float fxt;
	float fzo;
	float fxso;
	float fxgo;

	float fs;

	float dt;
	float dzt;

	float dxt;
	float dzo;
	float dxso;
	float dxgo;

	float ds;
	float ext;
	float ezt;

	float ezo;
	float exso;
	float exgo;

	float es;
	float s;
	float scal;

	float fxin;

	float **xsrf;
	float **zsrf;

	float v00;
	float v0;

	float dvz;
	float fmax;
	float angmax;
	float offmax;
	float rmax;
	float aperx;
	float sx;
	float gx;
	float nanglg;
	float dxi;
	float fxi;
	float ***dats,***ttab,***tb,***pb,***angb,**tsum,**tt,**szif,
                *nangl,*sz,**tbg,**pbg,**angbg,***sigb,**sigbg;
        Surface *srf;
	
	char  *infile="stdin",*outfile="stdout",*ttfile,*jpfile;
	FILE  *infp,*outfp,*ttfp,*jpfp,*hdrfp;

	int verbose=0;		/* verbose flag =1, show more info */

	/* hook up getpar to handle the parameters */
	initargs(argc, argv);
	requestdoc(1);

	/* open input and output files	*/
	if( !getparstring("infile",&infile)) {
		infp = stdin;
	} else  
		if ((infp=fopen(infile,"r"))==NULL)
			err("cannot open infile=%s\n",infile);
	if( !getparstring("outfile",&outfile)) {
		outfp = stdout;
	} else  {
		outfp = fopen(outfile,"w");
	}

	efseeko(infp,(off_t) 0,SEEK_CUR);
	efseek(outfp,(off_t) 0 ,SEEK_CUR);
	if( !getparstring("ttfile",&ttfile))
		err("must specify ttfile!\n");
	if ((ttfp=fopen(ttfile,"r"))==NULL)
		err("cannot open ttfile=%s\n",ttfile);
	if( !getparstring("jpfile",&jpfile)) {
		jpfp = stderr;
	} else  
		jpfp = fopen(jpfile,"w");

	/* get information from the first header */
	if (!fgettr(infp,&tr)) err("can't get first trace");
	nt = tr.ns;
	if (!getparfloat("dt",&dt)) dt = tr.dt/1000000.0; 
	if (dt<0.0000001) err("dt must be positive!\n");
	if (!getparfloat("ft",&ft)) ft = (float)tr.delrt/1000.; 
	
	/* get traveltime table parameters	*/
	if (!getparint("nxt",&nxt)) err("must specify nxt!\n");
	if (!getparfloat("fxt",&fxt)) err("must specify fxt!\n");
	if (!getparfloat("dxt",&dxt)) err("must specify dxt!\n");
	if (!getparint("nzt",&nzt)) err("must specify nzt!\n");
	if (!getparfloat("fzt",&fzt)) err("must specify fzt!\n");
	if (!getparfloat("dzt",&dzt)) err("must specify dzt!\n");
	if (!getparint("ns",&ns)) err("must specify ns!\n");
	if (!getparfloat("fs",&fs)) err("must specify fs!\n");
	if (!getparfloat("ds",&ds)) err("must specify ds!\n");
	ext = fxt+(nxt-1)*dxt;
	ezt = fzt+(nzt-1)*dzt;
	es = fs+(ns-1)*ds;

        if (!getparint("nxi",&nxi)) err("must specify nxi!\n");
        if (!getparfloat("fxi",&fxi)) err("must specify fxi!\n");
        if (!getparfloat("dxi",&dxi)) err("must specify dxi!\n");
        if (!getparint("sgn",&sgn)) err("must specify sgn!\n");

	/* optional parameters	*/
	if (!getparint("nxso",&nxso)) nxso = (nxt-1)*2+1;
	if (!getparfloat("fxso",&fxso)) fxso = fxt;
	if (!getparfloat("dxso",&dxso)) dxso = dxt*0.5;
        if (!getparint("nxgo",&nxgo)) nxgo = (nxt-1)*2+1;
        if (!getparfloat("fxgo",&fxgo)) fxgo = fxt;
        if (!getparfloat("dxgo",&dxgo)) dxgo = dxt*0.5;
        if (!getparint("antiali",&antiali)) antiali=1;
	if (!getparint("nzo",&nzo)) nzo = (nzt-1)*5+1;
	if (!getparfloat("fzo",&fzo)) fzo = fzt;
	if (!getparfloat("dzo",&dzo)) dzo = dzt*0.2;

	if (!getparint("verbose",&verbose))	verbose=0;

	fzo = fzo*sgn;
	dzo = dzo*sgn;

	exso = fxso+(nxso-1)*dxso;
	exgo = fxgo+(nxgo-1)*dxgo;
	ezo = fzo+(nzo-1)*dzo;
	if(fxt>fxso||fxt>fxgo||ext<exso||ext<exgo||fzt>fzo||ezt<ezo) 
		err(" migration output range is out of traveltime table!\n");

	if (!getparfloat("v0",&v0)) v0 = 1500;
	if (!getparfloat("dvz",&dvz)) dvz = 0;
	if (!getparfloat("angmax",&angmax)) angmax = 60.;
	if  (angmax<0.00001) err("angmax must be positive!\n");
	if (!getparfloat("aperx",&aperx)) aperx = 0.5*nxt*dxt;
	if (!getparfloat("offmax",&offmax)) offmax = 99999;
	if (!getparfloat("fmax",&fmax)) fmax = 0.25/dt;

	if (!getparint("ntr",&ntr))	ntr = 100000;
	if (!getparint("mtr",&mtr))	mtr = 100;
        checkpars();
        decodeSurfaces(sgn,&nsrf,&nxzsrf,&xsrf,&zsrf);
        makesurf(0.025*dxgo,nsrf,nxzsrf,xsrf,zsrf,&srf);

	fprintf(jpfp,"\n");
	fprintf(jpfp," Datuming parameters\n");
	fprintf(jpfp," ================\n");
	fprintf(jpfp," infile=%s \n",infile);
	fprintf(jpfp," outfile=%s \n",outfile);
	fprintf(jpfp," ttfile=%s \n",ttfile);
	fprintf(jpfp," \n");

        fprintf(jpfp," nzt=%d fzt=%g dzt=%g\n",nzt,fzt,dzt);
	fprintf(jpfp," nxt=%d fxt=%g dxt=%g\n",nxt,fxt,dxt);
 	fprintf(jpfp," ns=%d fs=%g ds=%g\n",ns,fs,ds);
	fprintf(jpfp," \n");
        fprintf(jpfp," nxi=%d fxi=%g dxi=%g sgn=%d\n",nxi,fxi,dxi,sgn);
        fprintf(jpfp," \n");
	fprintf(jpfp," nzo=%d fzo=%g dzo=%g\n",nzo,fzo,dzo);
	fprintf(jpfp," nxso=%d fxso=%g dxso=%g\n",nxso,fxso,dxso);
        fprintf(jpfp," nxgo=%d fxgo=%g dxgo=%g\n",nxgo,fxgo,dxgo);
	fprintf(jpfp," \n");

        /* compute Z coordinate for the surfaces  */
        szif = ealloc2float(2,nxi);
        sz = ealloc1float(ns);
        nangl = ealloc1float(nxi);
        zcoorSurfaces(fxi,dxi,nxi,fs,ds,ns,srf,szif,sz,nangl);

	/* compute reference traveltime and slowness  */
	rmax = MAX(es-fxt,ext-fs);
	rmax = MIN(rmax,0.5*offmax+aperx);
	nr = 2+(int)(rmax/dxgo);
        tb = ealloc3float(nzt,nr,ns);
        pb = ealloc3float(nzt,nr,ns);
        sigb = ealloc3float(nzt,nr,ns);
        angb = ealloc3float(nzt,nr,ns);
        for (is=0;is<ns;is++){
            v00 = v0 + sz[is]*dvz;
            timeb(nr,nzt,dxgo,dzt,fzt,sz[is],dvz,v00,tb[is],pb[is],sigb[is],
                angb[is]);
        }

	fprintf(jpfp," nt=%d ft=%g dt=%g \n",nt,ft,dt);
 	fprintf(jpfp," fmax=%g\n",fmax);
	fprintf(jpfp," v0=%g dvz=%g \n",v0,dvz);
 	fprintf(jpfp," aperx=%g offmax=%g angmax=%g\n",aperx,offmax,angmax);
 	fprintf(jpfp," ntr=%d mtr=%d \n",ntr,mtr);
	fprintf(jpfp," ================\n");
	fflush(jpfp);

	/* allocate space */
	dats = ealloc3float(nt,nxgo,nxso);
	ng = ealloc1int(nxso);
	ttab = ealloc3float(nzt,nxt,ns);
	tt = alloc2float(nzt,nxt);
	tsum = alloc2float(nzt,nxt);
        tbg = alloc2float(nzt,nr);
        pbg = alloc2float(nzt,nr);
        sigbg = alloc2float(nzt,nr);
        angbg = alloc2float(nzt,nr);

        memset((void *) dats[0][0],0,nxso*nxgo*nt*sizeof(float));

 	fprintf(jpfp," input traveltime tables \n");
                       
	/* compute traveltime residual	*/
	for(is=0; is<ns; ++is){
		nseek = (off_t) nxt*nzt*is;
		efseeko(ttfp,nseek*((off_t) sizeof(float)),SEEK_SET);
		fread(ttab[is][0],sizeof(float),nxt*nzt,ttfp);
		s = fs+is*ds;
		resit(nxt,fxt,dxt,nzt,nr,dxgo,tb[is],ttab[is],s);
	}

	fprintf(jpfp," start receiver datuming ... \n");
	fprintf(jpfp," \n");
	fflush(jpfp);
	
        mtmax = 2*dxgo*sin(angmax*PI/180.)/(v0*dt);
        if(mtmax<1) mtmax = 1;

	jtr = 1;
	s = fxso;
	ktr = 0;
	fxin = fxgo;
	for(is=0; is<nxso; ++is)
		ng[is]=0;

	/* Store headers in tmpfile while getting a count */
	hdrfp = etmpfile(); 

	do {

	    float as,res;
	    int is;
	    sx = tr.sx;
	    gx = tr.gx;

	    if(sx!=s) {
		s = sx;
                if(gx>fxgo) fxin = gx;
	    }
	
	    io = (int)((sx-fxso)/dxso);

	    if(MIN(sx,gx)>=fs && MAX(sx,gx)<=es && 
	       MAX(gx-sx,sx-gx)<=ABS(offmax) ){


                /* Number of traces in receiver's gathers  */
                ng[io]++;

                efwrite(&tr, 1, HDRBYTES, hdrfp);

		/*  Down or up ward continuation of the receivers  */
	    	as = (gx-fs)/ds;
	    	is = (int)as; 

		if(is==ns-1) is=ns-2;
		res = as-is;
		if(res<=0.01) res = 0.0;
		if(res>=0.99) res = 1.0;
		sum2(nxt,nzt,1-res,res,ttab[is],ttab[is+1],tsum);
                sum2(nr,nzt,1-res,res,pb[is],pb[is+1],pbg);
                sum2(nr,nzt,1-res,res,tb[is],tb[is+1],tbg);
                sum2(nr,nzt,1-res,res,sigb[is],sigb[is+1],sigbg);
                sum2(nr,nzt,1-res,res,angb[is],angb[is+1],angbg);

                as = (gx-fxi)/dxi;
                is = (int)as;
                if(is==nxi-1) is=nxi-2;
                res = as-is;
                if(res<=0.01) res = 0.0;
                if(res>=0.99) res = 1.0;
                nanglg = (1-res)*nangl[is] + res*nangl[is+1];

		dat2d(tr.data,nt,ft,dt,sx,gx,dats[io],aperx,nxgo,fxin,
		      dxgo,nzo,fzo,dzo,mtmax,gx,fmax,nxi,fxi,dxi,angmax,tbg,
                      pbg,angbg,nr,tsum,nzt,fzt,dzt,nxt,fxt,dxt,
                      antiali,sgn,szif,nanglg,sigbg,verbose);

	        ++ktr;
	        if((jtr-1)%mtr ==0 ){
			fprintf(jpfp," Datumed receiver %d\n",jtr);
			fflush(jpfp);
	    	}

	    }
	    ++jtr;
	} while (fgettr(infp,&tr) && jtr<=ntr);

	fprintf(jpfp," Datumed %d receivers in total\n",ktr);

	/* Seg-Y header					*/

	rewind(hdrfp);

	scal = 1/sqrt(2*PI)*dxgo;
	for(io=0; io<nxso; io++)  {
		for(ixo=0; ixo<ng[io]; ixo++) {

			efread(&tro, 1, HDRBYTES, hdrfp);
			memcpy((void *) tro.data,
			   (const void *) dats[io][ixo], nt*sizeof(float));

			for(it=0; it<nt; it++)
			   tro.data[it] *=scal;

			/* write out */
			fputtr(outfp,&tro); 

		}
	}

	fprintf(jpfp," \n");
	fprintf(jpfp," output done\n");
	fflush(jpfp);

	efclose(jpfp);
	efclose(outfp);
	efclose(hdrfp);

	free1int(ng);	    
	free2float(tsum);
	free2float(tt);
        free2float(tbg);
        free2float(pbg);
        free2float(sigbg);
        free2float(angbg);
        free3float(pb);
        free3float(tb);
        free3float(sigb);
        free3float(angb);
	free3float(ttab);
	free3float(dats);
	return(CWP_Exit());
}

/**********************************************************************
        Subroutines adapted from Dave Hale's modeling library

        decodeReflectors
        decodeReflector
        makeref                 Autor: Dave Hale, CSM, 09/17/91
**********************************************************************/

void decodeSurfaces(int sgn, int *nrPtr, int **nxzPtr, float ***xPtr,
		    float ***zPtr)
/*************************************************************************
decodeSurfaces - parse surfaces parameter string
**************************************************************************
Input :
sgn		Sign of the datuming process

Output:
nrPtr           pointer to nr an int specifying number of surfaces = 2
nxzPtr          pointer to nxz specifying number of (x,z) pairs defining the
                surfaces
xPtr            pointer to array[x][nr] of x values for each surface
zPtr            array[z][nr] of z values for each surface
**************************************************************************/
{
        int nr,*nxz,ir;
        float **x,**z;
        char t[4096],*s;

        /* count surfaces */
        nr = countparname("surf");
        if (nr==0) nr = 1;

        /* allocate space */
        nxz = ealloc1(nr,sizeof(int));
        x = ealloc1(nr,sizeof(float*));
        z = ealloc1(nr,sizeof(float*));

        /* get surfaces */
        for (ir=0; ir<nr; ++ir) {
                if (!getnparstring(ir+1,"surf",&s)) s = "0,0;99999,0";
                strcpy(t,s);
                if (!decodeSurface(t,sgn,&nxz[ir],&x[ir],&z[ir]))
                        err("Surface number %d specified "
                                "incorrectly!\n",ir+1);
        }

        /* set output parameters before returning */
        *nrPtr = nr;
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
}

/* parse one surface specification; return 1 if valid, 0 otherwise */

int decodeSurface (char *string, int sgn, int *nxzPtr, float **xPtr,
		   float **zPtr)
/**************************************************************************
decodeSurface - parse one surface specification
***************************************************************************
Input:
string          string representing surface
sgn		Sign of the datuming process

Output:
nxzPtr          pointer to number of x,z pairs
xPtr            array of x values for one surface
zPtr            array of z values for one surface
**************************************************************************/
{
        int nxz,ixz;
        float *x,*z;
        char *s,*t;

        s = string;
        s = strtok(s,",;\0");
        /* count x and z values, while splitting string into tokens */
        for (t=s,nxz=0; t!=NULL; ++nxz)
                t = strtok(NULL,",;\0");

        /* total number of values must be even */
        if (nxz%2) return 0;

        /* number of (x,z) pairs */
        nxz /= 2;

        /* 2 or more (x,z) pairs are required */
        if (nxz<2) return 0;

        /* allocate space */
        x = ealloc1(nxz,sizeof(float));
        z = ealloc1(nxz,sizeof(float));

        /* convert (x,z) values */
        for (ixz=0; ixz<nxz; ++ixz) {
                x[ixz] = atof(s);
                s += strlen(s)+1;
                z[ixz] = atof(s)*sgn;
                s += strlen(s)+1;
        }

        /* set output parameters before returning */
        *nxzPtr = nxz;
        *xPtr = x;
        *zPtr = z;
        return 1;
}

void makesurf(float dsmax,int nr,int *nu,float **xu,float **zu,
        Surface **r)
/*****************************************************************************
Make piecewise cubic surfaces
******************************************************************************
Input:
dsmax           maximum length of surface segment
nr              number of surfaces = 2
nu              array[nr] of numbers of (x,z) pairs; u = 0, 1, ..., nu[ir]
xu              array[nr][nu[ir]] of surface x coordinates x(u)
zu              array[nr][nu[ir]] of surface z coordinates z(u)

Output:
r               array[nr] of surfaces
******************************************************************************
Notes:
Space for the nu, xu, and zu arrays is freed by this function, since
they are only used to construct the surfaces.
*****************************************************************************/
{
        int ir,iu,nuu,iuu,ns,is;
        float x,z,xlast,zlast,dx,dz,duu,uu,ds,fs,rsx,rsz,rsxd,rszd,
                *u,*s,(*xud)[4],(*zud)[4],*us;
        SurfaceSegment *ss;
        Surface *rr;

        /* allocate space for surfaces */
        *r = rr = ealloc1(nr,sizeof(Surface));

        /* loop over surfaces */
        for (ir=0; ir<nr; ++ir) {

                /* compute cubic spline coefficients for uniformly sampled u */
                u = ealloc1float(nu[ir]);
                for (iu=0; iu<nu[ir]; ++iu)
                        u[iu] = iu;
                xud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,xu[ir],xud);
                zud = (float(*)[4])ealloc1float(4*nu[ir]);
                csplin(nu[ir],u,zu[ir],zud);

                /* finely sample x(u) and z(u) and compute length s(u) */
                nuu = 20*nu[ir];
                duu = (u[nu[ir]-1]-u[0])/(nuu-1);
                s = ealloc1float(nuu);
                s[0] = 0.0;
                xlast = xu[ir][0];
                zlast = zu[ir][0];
                for (iuu=1,uu=duu; iuu<nuu; ++iuu,uu+=duu) {
                        intcub(0,nu[ir],u,xud,1,&uu,&x);
                        intcub(0,nu[ir],u,zud,1,&uu,&z);
                        dx = x-xlast;
                        dz = z-zlast;
                        s[iuu] = s[iuu-1]+sqrt(dx*dx+dz*dz);
                        xlast = x;
                        zlast = z;
                }

                /* compute u(s) from s(u) */
                ns = 1+s[nuu-1]/dsmax;
                ds = s[nuu-1]/ns;
                fs = 0.5*ds;
                us = ealloc1float(ns);
                yxtoxy(nuu,duu,0.0,s,ns,ds,fs,0.0,(float)(nu[ir]-1),us);

                /* compute surface segments uniformly sampled in s */
                ss = ealloc1(ns,sizeof(SurfaceSegment));
                for (is=0; is<ns; ++is) {
                        intcub(0,nu[ir],u,xud,1,&us[is],&rsx);
                        intcub(0,nu[ir],u,zud,1,&us[is],&rsz);
                        intcub(1,nu[ir],u,xud,1,&us[is],&rsxd);
                        intcub(1,nu[ir],u,zud,1,&us[is],&rszd);
                        ss[is].x = rsx;
                        ss[is].z = rsz;
                        ss[is].c = rsxd/sqrt(rsxd*rsxd+rszd*rszd);
                        ss[is].s = -rszd/sqrt(rsxd*rsxd+rszd*rszd);
                }

                /* fill in surface structure */
                rr[ir].ns = ns;
                rr[ir].ds = ds;
                rr[ir].ss = ss;

             /* free workspace */
                free1float(us);
                free1float(s);
                free1float(u);
                free1float((float*)xud);
                free1float((float*)zud);

                /* free space replaced by surface segments */
                free1(xu[ir]);
                free1(zu[ir]);
        }

        /* free space replaced by surface segments */
        free1(nu);
        free1(xu);
        free1(zu);
}

/* Estimation of the Z coordinate of the surfaces  */
  void zcoorSurfaces(float fx,float dx,int nx,float fxs,float dxs,int nxs,
        Surface *srf,float **szif,float *sz,float *nangl)
/*****************************************************************************
From the cubic spline calculation, the Z coor. of the surface are estimated
******************************************************************************
Input:
fx              First x position in the surfaces
dx              x sampling interval
nx              number of output location points
fxs             x-coordinate of the first shot (travel-time tables)
dxs             x-coordinate increment of shots (travel-time tables)
nxs             number of shots (travel-time tables)
srf             surface structure

Output:
szif            z[] coordinates of the surfaces (output locations)
sz              z[] coordinates of the current surface (shot positions)
nangl           array of angles that the normal form with the vertical
*****************************************************************************/
{
        int   nss,ik,jx,nsi,nsf,ns,ixi,is,ix;
        float x,ax,ax0,az,temp;
        float **ssx,**ssz,*sss;
        SurfaceSegment *ss;

        nsi = srf[0].ns;
        nsf = srf[1].ns;
        ns = MAX(nsi,nsf);
        ssx = alloc2float(2,ns);
        ssz = alloc2float(2,ns);
        sss = alloc1float(nsi);

        /* loop over surface  */
        for(is=0; is<2; ++is){
            /* number of segments, segment length */
            nss = srf[is].ns;
            ss = srf[is].ss;
            for (jx=0; jx<nss ; jx++) {
                ssx[jx][is] = ss[jx].x;
                ssz[jx][is] = ss[jx].z;
                if(is==0) sss[jx] = ss[jx].s;
            }
            ixi = 0;
            for (ix=0; ix<nx; ++ix) {
                x = fx + ix*dx;
                for (ik=ixi; ik<nss-1; ++ik) {
                    if (ssx[ik][is]<=x && ssx[ik+1][is]>=x) {
                       az = ssz[ik][is] - ssz[ik+1][is];
                       ax0 = ssx[ik+1][is] - x;
                       ax = ssx[ik+1][is] - ssx[ik][is];
                       szif[ix][is] = ax0*az/ax+ssz[ik+1][is];
                       if (is==0){
                          az = sss[ik] - sss[ik+1];
                          temp = ax0*az/ax + sss[ik+1];
                          nangl[ix] = asin(temp);
                       }
                       ixi = ik;
                    }
                }
            }
            if(is==0){
                ixi = 0;
                for (ix=0; ix<nxs; ++ix) {
                    x = fxs + ix*dxs;
                    for (ik=ixi; ik<nss-1; ++ik) {
                        if (ssx[ik][is]<=x && ssx[ik+1][is]>=x) {
                            az = ssz[ik][is] - ssz[ik+1][is];
                            ax0 = ssx[ik+1][is] - x;
                            ax = ssx[ik+1][is] - ssx[ik][is];
                            sz[ix] = ax0*az/ax+ssz[ik+1][is];
                            ixi = ik;
                        }
                    }
                }
            }
        }
        for(is=0;is<2;is++){
            szif[0][is] = 2*szif[1][is] - szif[2][is];
            szif[nx-1][is] = 2*szif[nx-2][is] - szif[nx-3][is];
        }
        sz[0] = 2*sz[1] - sz[2];
        nangl[0] = 2*nangl[1] - nangl[2];
        sz[nxs-1] = 2*sz[nxs-2] - sz[nxs-3];
        nangl[nx-1] = 2*nangl[nx-2] - nangl[nx-3];

        free2float(ssx);
        free2float(ssz);
        free1float(sss);
}

/* residual traveltime calculation based  on reference   time	*/
  void resit(int nx,float fx,float dx,int nz,int nr,float dr,
		float **tb,float **t,float x0)
{
	int ix,iz,jr;
	float xi,ar,sr,sr0;

	for(ix=0; ix<nx; ++ix){
		xi = fx+ix*dx-x0;
		ar = abs(xi)/dr;
		jr = (int)ar;
		sr = ar-jr;
		sr0 = 1.0-sr;
		if(jr>nr-2) jr = nr-2;
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] -= sr0*tb[jr][iz]+sr*tb[jr+1][iz];
	}
} 

/* lateral interpolation	*/

/* sum of two tables	*/
  void sum2(int nx,int nz,float a1,float a2,float **t1,float **t2,float **t)
{
	int ix,iz;

	for(ix=0; ix<nx; ++ix) 
		for(iz=0; iz<nz; ++iz)
			t[ix][iz] = a1*t1[ix][iz]+a2*t2[ix][iz];
}
 
/* compute  reference traveltime and slowness	*/
      void timeb(int nr,int nz,float dr,float dz,float fz,float z0,float a,
	float v0,float **t,float **p,float **sig,float **ang)
{
	int  ir,iz;
	float r,z,v,rc,oa,temp,rou,zc;


	if( a==0.0) {
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz-z0;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				if(rou<dz) rou = dz;
				t[ir][iz] = rou/v0;
				p[ir][iz] = r/(rou*v0);
				sig[ir][iz] = v0*rou;
				ang[ir][iz] = asin(r/rou);
			}
	} else {
		oa = 1.0/a; 	zc = v0*oa;
		for(ir=0,r=0;ir<nr;++ir,r+=dr)
			for(iz=0,z=fz+zc-z0;iz<nz;++iz,z+=dz){
				rou = sqrt(r*r+z*z);
				v = v0+a*(z-zc+z0);
				if(ir==0){ 
					t[ir][iz] = log(v/v0)*oa;
					p[ir][iz] = 0.0;
					ang[ir][iz] = 0.0;
					sig[ir][iz] = 0.5*(z+0.1*dz-zc+z0)
						*(v0+v);
				} else {
					rc = (r*r+z*z-zc*zc)/(2.0*r);
					rou = sqrt(zc*zc+rc*rc);
					t[ir][iz] = log((v*(rou+rc))
						/(v0*(rou+rc-r)))*oa;
					p[ir][iz] = sqrt(rou*rou-rc*rc)
						/(rou*v0);
                                        temp = v0*p[ir][iz];
                                        if(temp>1.0) temp = 1.0;
                                        ang[ir][iz] = asin(temp);
					sig[ir][iz] = a*rou*r;
				}
			}
	}
}

void filt(float *trace,int nt,float dt,float fmax,int m,float *trf,
	     int sgn);

  void dat2d(float *trace,int nt,float ft,float dt,float sx,float gx,
	     float **dat,float aperx,int nx,float fx,float dx,
  	     float nz,float fz,float dz,int mtmax,float xm,
	     float fmax,int nxi,float fxi,float dxi,float angmax,
	     float **tb,float **pb,float **angb,int nr,float **tsum,
             int nzt,float fzt,float dzt,int nxt,float fxt,float dxt,
             int antiali,int sgn,float **szif,float nangl,
             float **sigb,int verbose)
/*****************************************************************************
Datum one trace 
******************************************************************************
Input:
*trace		one seismic trace 
nt		number of time samples in seismic trace
ft		first time sample of seismic trace
dt		time sampleing interval in seismic trace
sx,gx		lateral coordinates of source and geophone (ignored)
aperx		lateral aperture in the datuming process
nx,fx,dx,nz,fz,dz	dimension parameters of datuming region
mtmax		number of time samples in triangle filter
xm		datuming point ( source=sx or receiver=gx )
fmax		frequency-highcut for input trace	 
fxi             x-coordinate of the first surface location
dxi             horizontal spacingon surface
nxi             number of input surface locations
angmax		migration angle aperature from vertical 	 
(tb,pb,angb)	reference traveltime, lateral slowness, 
sigb		emergent angle and sigma
nr		number of lateral samples in reference quantities
tsum		sum of residual traveltimes from shot and receiver
nxt,fxt,dxt,nzt,fzt,dzt		dimension parameters of traveltime table
antiali         Anti-aliase filter flag
sgn             Sign of the datuming process(up or down ward)
szif            array[] of the z component of the reference surfaces
nangl		Angles that the normal unit vector of the
                current surface form with the vertical at source or geophone
                location

Output:
dat		Redatumed section
*****************************************************************************/
{
        int nxf,nxe,nxtf,nxte,it,ix,iz,iz0,izt0,nzp,jr,jz,jt,mt,jx,mr;
        float x,dis,rxz,ar,sr,sr0,z0,rdz,ampd,res0,am,am0,fxf,
	      ang,ax,ax0,pmin,odt=1.0/dt,pd,az,sz,sz0,
              at,td,res,temp,sig,sigma;
	float **tmt,**ampt,**ampti,*tm,*amp,*ampi,*tzt,*trf,*zpt;

	tmt = alloc2float(nzt,nxt);
	ampt = alloc2float(nzt,nxt);
	ampti = alloc2float(nzt,nxt);
	tm = alloc1float(nzt);
	tzt = alloc1float(nzt);
	amp = alloc1float(nzt);
	ampi = alloc1float(nzt);
	zpt = alloc1float(nxt);
	trf = alloc1float(nt+2*mtmax);

	z0 = (fz-fzt)/dzt + 0.0*sx + 0.0*gx;
	rdz = dz/dzt;
	pmin = 1.0/(2.0*dx*fmax);
	fxf = fxi + (nxi-2)*dxi;
	
	filt(trace,nt,dt,fmax,mtmax,trf,sgn);

	rxz = (angmax==90)?0.0:1.0/tan(angmax*PI/180.);
	nxtf = (xm-aperx-fxt)/dxt;
	if(nxtf<0) nxtf = 0;
	nxte = (xm+aperx-fxt)/dxt+1;
	if(nxte>=nxt) nxte = nxt-1;

	/* compute amplitudes and filter length	*/
	for(ix=nxtf; ix<=nxte; ++ix){
		x = fxt+ix*dxt;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = ((dis-dxt)*rxz-fzt)/dzt-1;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;

		ar = (xm>=x)?(xm-x)/dx:(x-xm)/dx;
		jr = (int)ar;
		if(jr>nr-2) jr = nr-2;
		sr = ar-jr;
		sr0 = 1.0-sr;
                sig = ((xm-x)<0)?1.0:-1.0;
		zpt[ix] = fzt+(nzt-1)*dzt;

		for(iz=izt0; iz<nzt; ++iz){
			sigma = sr0*sigb[jr][iz]+sr*sigb[jr+1][iz];
			ang = sr0*angb[jr][iz]+sr*angb[jr+1][iz];
                        ang = (sig*ang - nangl);
			ampd = cos(ang)/(cos(nangl)*sqrt(sigma));

	/* Filter of 90 degrees to the operator (Peels,1988)

           Peels, G. L., 1988, True amplitude wave field extrapolation
           with applications in  seismic shot record redatuming, PhD
           thesis, Delft University of Technology.              */

			if(ABS(ang)>=(PI/2)) {
				ampd=0.0;
			}

			if(ampd<0.0) ampd = -ampd;
			ampt[ix][iz] = ampd;

			pd = sr0*pb[jr][iz]+sr*pb[jr+1][iz];
			if(pd<0.0) pd = -pd;
			temp = pd*dx*odt;
			if(temp<1) temp = 1.0;
			if(temp>mtmax) temp = mtmax;
			ampti[ix][iz] = ampd/(temp*temp);
			tmt[ix][iz] = temp;
			if(pd<pmin && zpt[ix]>fzt+(nzt-1.1)*dzt) 
				zpt[ix] = fzt+iz*dzt;

		}
	}

	nxf = (xm-aperx-fx)/dx+0.5;
	if(nxf<0) nxf = 0;

	if((xm+aperx)>=fxf)
		nxe = (fxf-fx)/dx + 0.5;
	else
		nxe = (xm+aperx-fx)/dx + 0.5;

	if(nxe>=nx) nxe = nx-1;

	am = (fx-fxi)/dxi;

	mr = (int)am;
	am = am - mr;
	if(am<=0.01) am = 0.;
	if(am>=0.99) am = 1.0;
	am0 = 1.0-am;
	if(mr<0) mr = 0;
	if(mr+nxe>=nxi-1) 
		err("Topography definition is out of range!\n");
	
	/* interpolate amplitudes and filter length along lateral	*/
	for(ix=nxf; ix<=nxe; ++ix){
		x = fx+ix*dx;
		dis = (xm>=x)?xm-x:x-xm;
		izt0 = (dis*rxz-fzt)/dzt;
		if(izt0<0) izt0 = 0;
		if(izt0>=nzt) izt0 = nzt-1;
		iz0 = (dis*rxz-fz)/dz;
		if(iz0<0) iz0 = 0;
		if(iz0>=nz) iz0 = nz-1;

		ax = (x-fxt)/dxt;
		jx = (int)ax;
		ax = ax-jx;
		if(ax<=0.01) ax = 0.;
		if(ax>=0.99) ax = 1.0;
		ax0 = 1.0-ax;
		if(jx>nxte-1) jx = nxte-1;
		if(jx<nxtf) jx = nxtf;

		ar = (xm>=x)?(xm-x)/dx:(x-xm)/dx;
		jr = (int)ar;
		if(jr>nr-2) jr = nr-2;
		sr = ar-jr;
		sr0 = 1.0-sr;

		for(iz=izt0; iz<nzt; ++iz){
		    tzt[iz] = ax0*tsum[jx][iz]+ax*tsum[jx+1][iz]
				+sr0*tb[jr][iz]+sr*tb[jr+1][iz];

		    amp[iz] = ax0*ampt[jx][iz]+ax*ampt[jx+1][iz];
		    ampi[iz] = ax0*ampti[jx][iz]+ax*ampti[jx+1][iz];
		    tm[iz] = ax0*tmt[jx][iz]+ax*tmt[jx+1][iz];
		}

		nzp = (ax0*zpt[jx]+ax*zpt[jx+1]-fz)/dz+0.5;
		if(nzp<iz0) nzp = iz0;
                if(nzp>nz) nzp = nz;

                iz = (ABS(am0*szif[ix+mr][1]+am*szif[ix+mr+1][1])-fz)/dz;
                if(iz>=nz)
                        err("Datuming surface is out of output range!\n");
                az = z0+iz*rdz;
                jz = (int)az;
                if(jz>=nzt-1) jz = nzt-2;
                sz = az-jz;
                sz0 = 1.0-sz;
                td = sz0*tzt[jz]+sz*tzt[jz+1];
                at = (sgn*td-ft)*odt;
                if ((iz<nzp) && (antiali)) {
                /* interpolate along depth if operater aliasing   */
                        at = at + mtmax;
                        jt = (int)at;
                        ampd = sz0*ampi[jz]+sz*ampi[jz+1];
                        mt = (int)(0.5+sz0*tm[jz]+sz*tm[jz+1]);
                        res = ABS(at-jt);
                        res0 = 1.0-res;
                        for (it=0; it<nt; it++){
                          if(it+jt >= mtmax && jt < nt-it+mtmax-sgn){
                             temp = (res0*(-trf[it+jt-mt]+2.0*trf[it+jt]
                                 -trf[it+jt+mt])+res*(-trf[it+jt-mt+sgn]
                                 +2.0*trf[it+jt+sgn]-trf[it+jt+mt+sgn]))*ampd;
                             dat[ix][it] += temp;
                          }
                        }

                }
                /* interpolate along depth if not operater aliasing     */
                else{
                        jt = (int)at;
                        ampd = sz0*amp[jz]+sz*amp[jz+1];
                        res = ABS(at-jt);
                        res0 = 1.0-res;
                        for (it=0; it<nt; it++){
                          if(it+jt >= 0 && jt < nt-it-sgn){
                            temp=(res0*trace[it+jt]+res*trace[it+jt+sgn])*ampd;
                            dat[ix][it] += temp;
                          }
                        }

                }

	}

	free2float(ampt);
	free2float(ampti);
	free2float(tmt);
	free1float(amp);
	free1float(ampi);
	free1float(zpt);
	free1float(tm);
	free1float(tzt);
	free1float(trf);
}

void filt(float *trace,int nt,float dt,float fmax,int m,
	  float *trf,int sgn)
/*******************************************************************
  Low-pass filter, integration and phase shift for input data	 
  input: 
    trace(nt)	single seismic trace
    fmax	high cut frequency
    sgn		Sign of the Datuming process
  output:
    trace(nt) 	filtered and phase-shifted seismic trace 
    tracei(nt) 	filtered, integrated and phase-shifted seismic trace 
********************************************************************/
{
	static int nfft=0, itaper, nw, nwf;
	static float *taper, *amp, *ampi, dw;
	int  it,iw,itemp;
	float temp, ftaper, const2, *rt;
	complex *ct;

	fmax *= 2.0*PI;
	ftaper = 0.1*fmax;
	const2 = 0.5*sqrt(2.0);

	if(nfft==0) {
        	/* Set up FFT parameters */
        	nfft = npfaro(2*(nt+2*m), 4*(nt+2*m));
        	if (nfft >= SU_NFLTS || nfft >= 720720)
                	err("Padded nt=%d -- too big", nfft);

        	nw = nfft/2 + 1;
		dw = 2.0*PI/(nfft*dt);

		itaper = 0.5+ftaper/dw;
		taper = ealloc1float(2*itaper+1);
		for(iw=-itaper; iw<=itaper; ++iw){
			temp = (float)iw/(1.0+itaper); 
			taper[iw+itaper] = (1-temp)*(1-temp)*(temp+2)/4;
		}

		nwf = 0.5+fmax/dw;
		if(nwf>nw-itaper-1) nwf = nw-itaper-1;
		amp = ealloc1float(nwf+itaper+1);
		ampi = ealloc1float(nwf+itaper+1);
		amp[0] = ampi[0] = 0.;
		for(iw=1; iw<=nwf+itaper; ++iw){
			amp[iw] = sqrt(dw*iw)/nfft;
			ampi[iw] = 0.5/(1-cos(iw*dw*dt));
		}
	}

        /* Allocate fft arrays */
        rt   = ealloc1float(nfft);
        ct   = ealloc1complex(nw);

        memcpy(rt, trace, nt*FSIZE);
        memset((void *) (rt + nt), 0, (nfft-nt)*FSIZE); 
        pfarc(1, nfft, rt, ct);

	for(iw=nwf-itaper;iw<=nwf+itaper;++iw){
		itemp = iw-(nwf-itaper);
		ct[iw].r = taper[itemp]*ct[iw].r; 
		ct[iw].i = taper[itemp]*ct[iw].i; 
	}
	for(iw=nwf+itaper+1;iw<nw;++iw){
		ct[iw].r = 0.; 
		ct[iw].i = 0.; 
	}

	for(iw=0; iw<=nwf+itaper; ++iw){
		/* phase shifts PI/4 - Half derivative	*/
		temp = (ct[iw].r-sgn*ct[iw].i)*amp[iw]*const2;
		ct[iw].i = (sgn*ct[iw].r + ct[iw].i)*amp[iw]*const2;
		ct[iw].r = temp;
	}
              
        pfacr(-1, nfft, ct, rt);
		
        /* Load traces back in */
	for (it=0; it<nt; ++it) trace[it] = rt[it];

        /* Integrate traces   */
	for(iw=0; iw<=nwf+itaper; ++iw){
		ct[iw].i = ct[iw].i*ampi[iw];
		ct[iw].r = ct[iw].r*ampi[iw];
	}
        pfacr(-1, nfft, ct, rt);
        for (it=0; it<m; ++it)  trf[it] = rt[nfft-m+it];
        for (it=0; it<nt+m; ++it)  trf[it+m] = rt[it];

	free1float(rt);
	free1complex(ct);
}
